function [T_p, T_l_p, T_u_p] = M_p( ln_m, vec)

    %%
    %% First derivative
    %%

    n=length(ln_m);

    ind_im = 1+1 + (n+1).*[0:n-2];
    ind_ip   = 1-1 + (n+1).*[1:n-1];
    ind_i = 1   + (n+1).*[0:n-1];

    T_u_p=sparse(n,n);
    T_l_p=sparse(n,n);

    dinv=1./(ln_m(2:end)-ln_m(1:end-1));

    el_i   = ones(n,1).*(1./2.*[nan(1);dinv] ...
        -1./2.*[dinv;nan(1)]);

    el_i_u   = ones(n,1).*(-[dinv;nan(1)]);

    el_i_l   = ones(n,1).*([nan(1);dinv]);

    el_i_u(1)=-dinv(1);
    el_i_u(end)=dinv(end);

    el_i_l(1)=-dinv(1);
    el_i_l(end)=dinv(end);

    el_i(end)=dinv(end);
    el_i(1)=-dinv(1);

    el_im =   -ones(n-1,1).*(dinv);
    el_im(end)=-dinv(end);
    el_ip   =   ones(n-1,1).*(dinv);
    el_ip(1) =   dinv(1);

    el_im_u =   zeros(n-1,1);
    el_im_u(end)=-dinv(end);
    el_ip_l   =   zeros(n-1,1);
    el_ip_l(1) =   dinv(1);

    T_u_p(ind_im) =T_u_p(ind_im)'+ el_im_u.*vec(1:end-1);
    T_u_p(ind_ip) = T_u_p(ind_ip)'+ el_ip.*vec(2:end);
    T_u_p(ind_i) = T_u_p(ind_i)'+  el_i_u.*vec;

    T_l_p(ind_im) = T_l_p(ind_im)'+ el_im.*vec(1:end-1);
    T_l_p(ind_ip) = T_l_p(ind_ip)'+ el_ip_l.*vec(2:end);
    T_l_p(ind_i)  = T_l_p(ind_i)'+  el_i_l.*vec;

    T_p = (T_l_p+T_u_p)./2;

end
